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Abstract 

An one - parameter dynamical system is associated to the mathe- 
matical problem governing the membrane excitability of a ventricular 
cardiomyocyte, according to the Luo-Rudy I model. Limit cycles are 
described by the solutions of an extended system. A finite element 
method time approximation (FEM) is used in order to formulate the 
approximate problem. Starting from a Hopf bifurcation point, ap- 
proximate limit cycles are obtained, step by step, using an arc-length- 
continuation method and Newton's method. Some numerical results 
are presented. 

Key words: limit cycle, finite element method time approximation, 
Luo-Rudy I model, arc-length-continuation method, Newton's method. 
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1 Introduction 



The well-known Hodgkin-Huxley model of the squid giant axon ([16j) 
represented a huge leap forward compared to earlier models of excitable sys- 
tems built from abstract sets of equations or from electrical circuits including 
non-linear components, e.g. [33] , The pioneering work of the group of Denis 
Noble made the transition from neuronal excitability models, characterized 
by Na-|- and K-|- conductances with fast gating kinetics, to cardiomyocyte 
electrophysiology models, a field expanding steadily for over five decades 
(|23)). Nowadays, complex models accurately reproducing transmembrane 
voltage changes as well as ion concentration dynamics between various sub- 
cellular compartments and buffering systems are incorporated into detailed 
anatomical models of the entire heart (|24j). The Luo-Rudy I model of iso- 
lated guinea pig ventricular cardiomyocyte (|21j) was developed in the early 
1990s starting from the Beeler- Renter model ([I]). It includes more recent 
experimental data related to gating and permeation properties of several 
types of ion channels, obtained in the late 1980s with the advent of the 
patch-clamp technique ([22]). The model comprises only three time and 
voltage-dependent ion currents (fast sodium current, slow inward current, 
time-dependent potassium current) plus three background currents (time- 
independent and plateau potassium current, background current), their dy- 
namics being described by Hodgkin-Huxley type equations. This apparent 
simplicity, compared to more recent multicompartment models, renders it 
adequate for mathematical analysis using methods of linear stability and 
bifurcation theory. 

Nowadays, there exist numerous software packages for the numerical 
study of finite - dimensional dynamical systems, for example MATCONT, 



CL_MATCONT, CL.MATCONTM ([7], [E]), AUTO [8]. In [19], [8], [T], 



[15] , the periodic boundary value problems used to locate limit cycles are ap- 
proximated using orthogonal collocation method. Finite differences method 
is also considered. In this paper, limit cycles are obtained for the dynamical 
system associated to the Luo-Rudy I model by using finite element method 
time approximation (FEM). 

2 Luo-Rudy I model 

The mathematical problem governing the membrane excitability of a 
ventricular cardiomyocyte, according to the Luo-Rudy I model ([21]), is a 
Cauchy problem 



u{0) = no , 

for the system of first order ordinary differential equations 



(1) 




(2) 
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where u = {ui,...,us) = {V, [Ca\i, h, j, m, d, /, X), rj = (771, . . . , r^is) 
= {1st, Cm, QNa, 9si, QKp, Qb, [Na]^, [Na]i, [K]o, [K]i, PRNaK, Eb, T), 
M = R^, T -.R^^ X M ^ M, T = {Ti,...,Ts), 

1 „ 

Ti{r],u) = [1st +mU3U4U^{ui - ENaim, 118,^13)) 

m 

+riiUQU7{ui - ci + C2 111^2) 

+9K{mo)Xi{ui){ui - EK{m,V8,V9,mo,Vu,Vn))us 
+9Ki{vw)Kloc{r]9,Vw,Vi3,ui){ui - mo, ms)) 

+r]5Kp{ui){ui - Ekp{v9, mo,Vri)) + Vei^i - V12)] , 

J^2{'n,u) = -C3r]4UQU7{ui -Ci +C2lnU2) +C4(C5 - ^2) , 

Ti{7], u) = ai{ui) - {ai{ui) + I3i,{ui))ui, , i = 3,...,8. 

For the definition of variables V, [Ca]i, h, j, in, d, f, X, parameters 1st, 
Cm, 9Na, 9 si, 9Kp, gb, [Na]o, [Na]i, [K]q, [K]i, PRNaK, Eb, T, constants 
ci, . . . , C5, functions qk, ENa, Ek, Eki, Exp, i^loo, Xi, Kp, a^, (3i>, default 
values of parameters and initial values of variables in the Luo-Rudy I model, 
the reader is referred to [21]. The reader is also referred to [20j for the 
continuity of the model, and to [1] for the treatment of the vector field J- 
singularities. J- is of class C^ on the domain of interest. 

3 The one - parameter dynamical system associ- 
ated to the Luo - Rudy I model 

We performed the study of the dynamical system associated with the 
Cauchy problem ([I]), ([2]) by considering only the parameter r/i = 1st and 
fixing the rest of parameters. Denote A = ?7i = 1st and r/* the vector of 
the fixed values of r/2, . . . , r/13. Let F : M x M — )• M, -F(A, u) = F[\, ri^,,u), 
F = {Fi,...,Fs). 

Consider the dynamical system associated with the Cauchy problem ([1]), 
([3]), where 

du . 

^ = F(A.„). (3) 

The equilibrium points of this problem are solutions of the equation 

F{X,u)=0. (4) 

The existence of the solutions and the number were established by graphical 
representation in [3], for the domain of interest. The equilibrium curve 
(the bifurcation diagram) was obtained in [3] , via an arc- length-continuation 
method ([13j) and Newton's method ([12j), starting from a solution obtained 
by solving a nonlinear least-squares problem ( [13] ) for a value of A for which 
the system has one solution. In [1] , the results are obtained by reducing (H]) 
to a system of two equations in (ui, M2) = (^1 [C'«]j)- Here, we used directly 
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4 Extended system method for limit cycles 

The extended system in (A, T, u) 

^-TFiX,u) = 0, 

u(0)-n(l) = 0, (5) 
j <uit),^>dt = 



is introduced, in [12], [8], [7], to locate limit cycles of a general problem ([T|), 
([3]). T is the unknown period of the cycle, w is a component of a known 
reference solution {X,T,w) of ([5]). The system ([5]) becomes determined in a 
continuation process. 

8 ri 

In our case, < u,v > = ^i^i |[n|| = < / ^ n? for ti, u € M*. 

1=1 y i=i 

In order to approximate and solve ^ by finite element method time 
approximation (FEM), let us obtain the weak form of ([5]) in the sequel. 
Let 

X = {x € L2(0, 1; M^); ^ G ^^(0, 1; M*^), 

x = {xi,..., xs), Xi(0) = i = 1, . . . , 8} . 

The weak form of ([5|) is the problem in (A, T, n) € M x M x X 

1 1 

/ ^^dT + T J F,{X,u{T))v{T)dT = 0, 



yveV,i = l,...,8, (6) 

1 



< u(t), — ^ >dt = 0. 
at 



5 Arc-length-continuation method for ([6]) 

Following the usual practice ([IT], [IB], [7], [8], [12], [l3], [H], [l5], [l9], 
[25] , ^27j , ^22] , [2U] ) , we also use an arc-length-continuation method in order 
to formulate an algorithm to solve Q approximatively. 

Glowinski ([E], following H.B.Keller [TT], [12]) and Doedel (0, where 
also Keller's name is cited) chose a continuation equation written in our case 
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as 



, du(t) ,,9 , ,dT .r, ,(iA, o 
' ^ '''^dt + {—f + {—f 



ds 



ds 



■ ds ' 



1, 



(7) 



where s is the curvihnear abscissa. 

Let {}^,v^) be a Hopf bifurcation point, a pair of purely imag- 

inary eigenvalues of of the Jacobian matrix DuF{X^ ,u^), and a nonzero 
complex vector = gr + Wi- {X^,u^) is located on the equilibrium curve 
during a continuation procedure using some test functions ([IH], [H], [7]). 

system ([27], [SH], [29]) 



X 11^^ X 



is the solution of the extended 



F{X,u) 

DuF{\, u)gr + Pgi 
DuF{\, u)gi - (3gr 

gr,k — 1 

9i,fc 



0, 



(8) 



where A: is a fixed index of gr and of (7j, 1 < A; < 8. 

To solve ([6]), the extended system formed by ^ and d?]), parametrized 
by s, was considered. Let As be an arc-length step and A" = u(AAs), 
rpn ^ T(nAs), vP' = u{nAs). We have the algorithm (following the cases 
from [E], [S], [^, r29j): 

1. take the Hopf bifurcation point {X^,u^) and = 27r//3'^; retain (^J!, 



5°; 



2. for n = 0, (A\r\ni) G 



X X is obtained ([8], [29]) by (fTHD . 

#i(t) 



dt 



dt = , 



(9) 



and 



(10) 



where 

(j){t) = sin(27rt)c/° + cos(27rt)5r° , 
using Newton's method with the initial iteration 

{uy{t) = u^ + Ascj){t), {Ty = T\ 



{X')" = A" 



(11) 
(12) 
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3. for n > 1, assuming that (A"-\ r"-\ li""^), (A", T", u") are known, 
(A"+\ T"+i, e M X R X X is obtained by (USD, (HH), and where 



1 

dr + T^+^ I Fi{X^+\u^+\r))v{T)dT = 0, (13) 
V^; G i = 1,...,8, 



1 8 

lY^ur\t)^dt=o, (14) 



^«+^(t) - n?(t))<^^l^^ dt (15) 

i=l 

T-in rpn—l \n X""! 

- T'')- + (A"+i - A")^-— ^ = As , 

As As 

using Newton's method with the initial iteration 

((A"+i)0, (r"+i)0, (n"+i)°) = (A", T", n"). . (16) 

6 Newton's method for the steps of the algorithm 
from the end of section [5] 

In ([IS]) (n > 1), let us denote A* = A", T* = T", u* = u"", A** = 
A!!^, T- = = We write m, ©, (ffOl) (the 

iteration n = 0) in the same general form as (fT3]l . (fT^ . (fTSj) . So denote u* 
= ti** = (p and consider A* = A°, T* = T^, A** = 0, T** = in ([l5]) 
and consider u* = = ip in ()14p . 

Each step of the algorithm at the end of section [5l given (A*, T*, u*), 
(A**, T**, u**), calculates (A'^+i, r"+i, u"+i) G M x M x X, n > 0, by 

1 g 

lE<^'(')^dt = 0, (17) 




and 



8 

[^{urHt)-uut))ur{t)dt (18) 



+ (^"+1 _ rp*^'J-** ^ (_xn+i _ _x*)A** = As . 
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Newton's method applied 1^, ^ and ([T8]), for n > 0, leads to: let 
((A^)°, (Ti)0, given by ([H]), be an initial iteration (m = 0) if n = 0; 

let ((A"+^)°, (r"+^)°, given by ([IS]), be an initial iteration (m = 0) 

if n > 1; calculate (A"^-'^, T""'"^, u""'"^) as the solution of the algorithm: for 
m > 0, ((A"+i)™+i, (r"+i)"^+i, = (A"^+\ r"^+\ n™+i) G R x 

M X X is obtained by 

1 1 

J ^ dT + T^+' I F,{\^,U^{T))v{T)dT, 



1 

J DF,(A™,n'"(T))(A"+\tx™+i(r))t;(r)(ir (19) 





1 



_|_ ji m 




1 



T"" J DFi{y^, n"^(T))(A™, n™(T))i;(T)dr 



G i = 1, 



L, . . . , W , 



E-r'(*)^^* = 0' (20) 



/ dt + T'"+iT** + A"+U** (21) 

\ 8 

/ E u* {t)u** (t) dt + T*T** + A* A** + As . 







i=l 







i=l 



7 Approximation of problem (fT9l). (1201). ( 1211 ) by 



finite element method time approximation 

In order to perform this approximation, let us divide the interval [0, 1] 
in iV + 1 subintervals K = Kj = [tj,tj+i], < j < N, where = to < ti < 
. . . < tAT+i = 1. The sets K represent a triangulation Th of [0, 1]. 

Let us approximate the spaces V and X by the spaces 

Vh = {v: [0, 1] ^ M; t; G C[0, 1], ^;(0) = t;(l), v\k G Pfc(^), VA^ G Th} , 
Xh = {x: [0,1] ^M^; X = (xi,...,X8), Xi G i = 1, . . . , 8} , 

respectively, where Pk{K) is the space of polynomials in t of degree less than 
or equal to k defined on JC, A: > 2. 
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Let k = 2. An element K & Th has three nodal points. To obtain a 
function G reduces to obtain a function G Vh- In order to obtain 
a function Vh S Vh, we use a basis of functions of Vh- Let Jk = {1, 2, 3} be 
the local numeration for the nodes of K, where 1,3 correspond to tj,tj+i 
respectively and 2 corresponds to a node between tj and tj+i- Let {tpiji £ 
Jk} be the local quadratic basis of functions on K corresponding to the 
local nodes. Let J = {1, . . . , 2N + 1} be the global numeration for the nodes 
of [0, 1]. The two numerations are related by a matrix L whose elements are 
the elements j G J. Its rows are indexed by the elements K €z Th (by the 
number of the element K m a certain fixed numeration with elements from 
the set {1, ... , A^}) and its columns, by the local numeration i £ Jk, that is 
j = L{K,i). A function Vh € Vh is defined by its values Vj from the nodes 
j e J, 

Mt)= Yl E ^jMt), (22) 

K& Th ie JK,j=L{K,i) 

and a function Uh G Xh is defined by its values Uj from the nodes j € J, 

Mt)= Yl Yl u,Mt)- (23) 

A'e Th i€ JK,j=L{K,i) 
So, an unknown function = {iuh)i, ■ ■ ■, (iih)8) is reduced to the unknowns 

Uj, Uj = {{Uj)l, . . ., {Uj)s), j G J. 

In ([IS]), dZOD, ([21]), approximate (A'^+\ r™+i, u™+i) E M x M x X by 
(A™+\ r™+i, u'^^^) eWxRx Xh. Taking = -u,,, n/, given by ([231), 

and V = Tpi, for all ^ € J_ft:, for all X G Th, we obtain the discrete variant of 
problem (fTUj) . ([201) . ([21]) as the following problem in (A, T, ui, . . ., U2n+i) 
G M X M X written suitable for the assembly process, 



Ke Th ie JK,j=L{K,i) ^ 

+ r-+i / Fi(A'",<(^))^^^^ 



+ Y < (^«i^n(A'",<(^))V'^(T)) dr > 

«e -/k, j=L{K,i) ^ 

+A™+i y (DAi^n(A™,^z;r(^))) V'^ } (24) 

K 

^ { j{DM^^,ut{r))u^{r))^ar 
Kd Th K 

+ / (I)Ai^n(A"^, n;r(^))A") V'^ } , n = 1, . . . , 8 , 



K 

m+1 



(25) 
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E E <nY^\lMr)^dr>=0, (26) 

'e Th i£ Jk, j=L{K,i) ^ 



Th ie Jk, j=L{K,i) 
KG Th ie Jk, j=L{K,i) ^ 



(27) 



Ke ThK 

for all ^ e Ji^, for all K € 7^ 



E / < > dT + r*T** + A*A** + As, 



8 Numerical results 

Based on [30j and on the computer programs for [2] and [3], relations 
(pi]) . (f25]) . (p6]) . ([27]) and the algorithm at the end of section [5] furnished the 
numerical results of this section. 

Let (A'^, u^) be the Hopf bifurcation point located during the construction 
of the equilibrium curve by a continuation procedure in [S]. The solution 
(A°,/3°,u°,5°,c/°) of (IHD, calculated in [5], is A° = -1.0140472901, /3° = 
0.0162886062, = (-24.3132508542, 0.0034641214, 0.0, 0.0, 0.9176777444, 
0.5025242162, 0.4920204612, 0.5071561613), 5° = (1.0, 0.0000468233, 0.0, 
0.0, 0.0093195354, 0.0198748652, -0.0072420216, 0.0001706577), = (0.0, 
0.0000029062, 0.0, 0.0, -0.0000171311, -0.0118192370, 0.0136415789, -0.0017802907). 
(The eigenvalues of the Jacobian matrix DuF{\^ ,u^), calculated by the QR 
algorithm, are ±0.0162886062 i, -8.8611865338, -0.1026761869, -0.0647560667, 
-0.0024565181, -1.7398266947, -0.2049715178). These data are considered 
in the step 1 of the algorithm at the end of section [SJ 

We took C„^ = 1, gNa = 23, gsi = 0.09, = 0.282, qki = 0.6047, 
gxp = 0.0183, Gb = 0.03921, [Na\o = 140, [Na]i = 18, [K]o = 5.4, [K]i = 
145, PRNaK = 0.01833, Eb = -59.87, T = 310. 

In order to solve ^ numerically by the algorithm at the end of section 
[5] and by ([MD, (ES]), we performed calculations using As = 1.0 

and 500 iterations in the continuation process. Integrals J /(t) dr were 

K 

calculated using Gauss integration formula with three integration points. 

Figure [1] and [2] present some results obtained using 20 elements K (41 
nodes) (A^ = 20, J = {1, . . . , 41} in section[7]). The curves of the projections 
of the limit cycles, on the planes indicates in figure, are plots generated 
from values calculated in the nodes, corresponding to a fixed value of the 
parameter. 

Two projections of some limit cycles and of a part of the equilibrium 
curve (marked by "A") are presented in figure [TJ The Hopf bifurcation 
point is marked by 
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In figure [21 there are represented the projections of the plots of two 
Umit cycles calculated for 1st = —1.2000465026 (iteration 148) and for Igt = 
-1.2000183729 (iteration 248, marked by "x" in figure). 

The results obtained are relevant from a biological point of view, point- 
ing to unstable electrical behavior of the modeled system in certain condi- 
tions, translated into oscillatory regimes such as early afterdepolarizations 
([32j) or self-sustained oscillations ([3]), which may in turn synchronize, re- 
sulting in life-threatening arrhythmias: premature ventricular complexes or 
torsades-de-pointes, degenerating in rapid polymorphic ventricular tacycar- 
dia or fibrillation (|26)). 
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